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We present efficient circuits that can be used for the phase space tomography of quantum states. 
The circuits evaluate individual values or selected averages of the Wigner, Kirkwood and Husimi 
distributions. These quantum gate arrays can be programmed by initializing appropriate computa- 
tional states. The Husimi circuit relies on a subroutine that is also interesting in its own right: the 
efficient preparation of a coherent state, which is the ground state of the Harper Hamiltonian. 
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I. INTRODUCTION 

Phase space distributions have been used as represen- 
tation tools for quantum mechanical operators since the 
early days of quantum mechanics. They provide the ideal 
link to explore and understand the transition to classi- 
cal mechanics and to display in phase space quantum 
effects. Their properties are very well known [HQ when 
the phase space is R 2 . For systems with a finite dimen- 
sional space of states the distributions become discrete, 
i.e. they are defined over a finite lattice 0,3- Discrete 
phase space distributions have been used in the context 
of studies of quantum maps on bounded phase space 0,0 
and they have also recently proposed as a useful tool for 
studies related to quantum information and computation 
0, IE 13 The simplest way to characterize them, for a 
Hilbert space of dim N is by using a complete basis of 
iV 2 operators {A a ;a = 0, ..N 2 — 1}, in terms of which 
the distribution is given as TrbA^p]. The properties and 
classical features that these distributions display depend 
of course on the operator basis {A a }. In this sense phase 
space distributions are nothing but the coefficients of the 
expansion of the state p in the basis A Q . The deter- 
mination of the value of Tr[A Q( o] for every a is, thus, a 
particular form of quantum state tomography (see [lOj 
and references therein). 

In a recent paper it was shown how to efficiently 
measure the discrete Wigner function at any phase space 
point. The basis of the method is the use of the so-called 
'scattering circuit' to efficiently determine the value of 
the quantity Tr[pyl] provided that the operation A can be 
implemented in a controlled way. Thus, if the complete 
basis A a consists of unitary operators then the scatter- 
ing circuit can be used to measure individual values of 
the distribution. The disadvantage is of this approach 
is that, as a enters as a classical parameter, a new gate 
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array has to be applied for each a. In this paper we will 
extend the results presented in 0O] in two ways. First, 
we will show how to efficiently measure other phase space 
distribution functions (Husimi, Kirkwood). Second, we 
will show how to do this by using quantum circuits with 
a fixed architecture, which is independent of the phase 
space point a. These circuits belong to the class of pro- 
grammable quantum devices, whose action is controlled 
by quantum software, that have been under investigation 
recently [l2l InH ] . In this paper we will show how to build 
efficient programmable circuits to measure three phase 
space distributions: Wigner, Kirkwood and Husimi. It is 
also worth mentioning here that the quantum circuits we 
developed use a subroutine which is interesting in its own 
right and could be useful for other applications. In fact, 
in this paper we present a method to efficiently prepare 
coherent sates (which are rigorously defined below, but 
can be roughly characterized as approximately Gaussian 
wave packets obeying periodic boundary conditions). 

The paper is organized as follows in Section [H] we 
present the circuit that enables the programmable mea- 
surement of the discrete Wigner function. We also show 
that it can be useful to compute averages of this function 
over various phase space domains (this extends and com- 
pletes results presented in 01). In Section llTTI we present 
a simple programmable circuit that evaluates the Kirk- 
wood distribution at any phase space point. In Section lTvl 
we present the quantum gate array that efficiently eval- 
uates the discrete Husimi distribution. This gate uses 
coherent states as inputs. The algorithm to efficiently 
prepare those states is presented in Section Finally, 
we present some conclusions in Section IVll 



II. PROGRAMMABLE TOMOGRAPHY OF 
THE DISCRETE WIGNER FUNCTION 

The discrete Wigner function 0, 0| in a Hilbert space 
of dimension N is defined in terms of the basis of phase 
point unitary operators: 

A(q,p) = U q RV~ p e l ^ pq , (1) 
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W(q,p) = — Tr[A(q,p)p] 



(2) 



where q,p are integer labels spanning a grid of size 2N x 
2N. U and V are respectively the translation operators 
in the \q) and \p) basis (U\q) = \q+l), U\p) = e~ i! fr p \p), 

V\p) = \p + l),V\q) = e% ~^ q \<l))i which are related by the 
discrete Fourier transform. R is the reflection operator 
(R\n) = \N — n)). Only a,n N x N sub-grid is needed 
for the complete tomography of the state (but the larger 
grid is required to define a Wigner function with all the 
desired properties @). 

The programmable circuit implementing the measure- 
ment of the discrete Wigner function is shown in Figure 
^ This was introduced in [b| as a particular case of a 
programmable circuit evaluating the expectation value of 
an arbitrary operator. It is a variation of the so-called 
scattering circuit |ll| where an ancillary qubit acts as 
a probe for a more complex system with which it in- 
teracts by means of controlled operations. The circuit 
shown in Figure Q has several registers: The first regis- 
ter is an ancillary qubit (the probe) initially prepared in 
the state |0) which is an eigenstate of a z with eigenvalue 
+1. The following two registers act as program regis- 
ters and should be prepared in the state \q)\p)- The last 
register stores the state of the system of interest p. The 
program state contains the information about the binary 
expansion of the coordinates of the point in the phase 
space where we wish evaluate the Wigner function. As 
seen in the circuit, the role of the program states is to 
control the application of displacement operators on the 
system register. Here, and in what follows, we use the 
convention that for any operator O, "controlled-O" op- 
erators act as: (ctrl-0)|n) |*) = \n)O n \^). In particular, 
in Figure ^ an operator such as "control- V~2n" acts as 
(ctrl-VW)j<z)b) = \q)V 2 q N \p) = e^™\q)\p) (note that a 
subscript in any operator indicates the dimensionality of 
the space in which it acts) . It is straightforward to show 
that the final polarization of the ancillary qubit turns out 
to be: 

(a z ) = 2NTi[A(q,p)p] = 2NW(q,p). (3) 



A. Measuring the sum of the discrete Wigner 
function over domains in phase space 

One of the defining properties of the Wigner function 
is the fact that adding its values over lines in phase space 
one always obtains the probability to measure an observ- 
able. It is interesting to notice that the circuit shown 
in Figure ^ can be programmed to directly evaluate the 
average of the Wigner function along any line in phase 
space. More generally, the state of the program register 
can be used to define the phase space domain over which 
the Wigner function is averaged. 
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FIG. 1: Programmable circuit to evaluate the discrete Wigner 
function from the polarization of the first qubit. The second 
and the third registers store the information about the the 
phase space point where the distribution is being measured. 
All states are initialized in the computational (coordinate) ba- 
sis. All "ctrl-O" operations act as (ctrl-0)|n)|*) = |n)O n |*). 
A subscript in an operator denotes the dimension of the space 
in which it acts. In all figures we adopt the following graphic 
convention: thin cables denote a single qubit, and wide cables 
denote systems of more than one qubit. 



Let us consider first the case of lines. The quantity 
in which we are interested is Yliq p)eL W(q,p), the sum 
of the values of the Wigner function along the line L. 
It is easy to see that for the program state \^)p = 
Ep=o 1 l n s) ® Ip)/V^> the final polarization turns out 
to be (a z ) = J2 P W(ri3 7 p). As this type of program state 
can be efficiently constructed, the example shows that 
it is possible to estimate the sum of the values of the 
Wigner function along vertical and horizontal lines. In a 
recent work 0] we showed that this is a special case of 
a more general result that establishes the possibility to 
program the measurement of the expectation value of any 
operator. Following the same idea, consider the program 
state 



{q',p')eL 



(4) 



where L : n\q + ri2p = ri3, < rii < 2N — 1 and K is 
a normalization constant (the square root of the number 
of points in the line L). Then, the final polarization is: 



= i E W(q', P '). 



E 

(q'.p')eL 



(•5) 



This is precisely the quantity we are interested in. How- 
ever, the above program states may be difficult to prepare 
(they are, in general, highly entangled states). To avoid 
using program states which may be difficult to prepare we 
have developed an alternative method. This was briefly 
described in |13| . For completeness, we present it here in 
more detail. 

For our method it is convenient to employ the fact 
that certain unitary operators induce a purely classical 
transformation of the Wigner function (this means that 
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the Wigner function is simply transported by a canoni- 
cal, area preserving, flow). This is the case for unitary 
operators that quantize linear canonical transformations 
on the torus (the so called cat maps [t| ) . We can use this 
fact as follows: First, we can prepare a simple program 
state that would produce the measurement of the Wigner 
function along vertical or horizontal lines. Then, we can 
obtain the corresponding measurement along tilted lines 
by applying the appropriate unitary cat map to the ini- 
tial state. The two-parameter family of cat operators 
that we will use is given by: 

Ucat = V b TV ai (6) 

where a and b are integers, and the operators V a and T 
are diagonal in the position and momentum basis respec- 
tively, 

V a \n) = exp(-i27m 2 (l - a)/2N)\n) 
T\k) = cx-p(~i2nk 2 /2N)\k). (7) 

The classical equations of motion corresponding to this 
system are 

q = bq'+p' p = (ab — l)q + ap . (8) 

As the Wigner function evolves classically, when the 
phase space points are related as above, we can write 
W(q,p, t+1) = W(q',p',t). The transformation JSJ) maps 
vertical lines into tilted lines according to the values of 
the parameters a and b. With this in mind we can try to 
find the linear transformation that maps a line L, whose 
program state is difficult to prepare, into a line L\ whose 
program state is easy to prepare. If we achieve this, we 
can compute the average Wigner function along L by 
using the fact that 

£ W p (q,p)= £ W Ucatpul Jq',p'). 
(q,p)eL ( q >,p>)eL> 

The circuit to implement this procedure is shown in Fig- 
ure |21 In order to implement the U ca t evolution, we only 
need to know how to apply the unitary operator T and 
its powers. As it was shown in this can be done ef- 
ficiently. Taking into account that the parameters defin- 
ing the evolution U ca t depend on the line, this network is 
not completely programmable (its architecture depends 
on each line), but we can easily prove that the device can 
transformed into a fully programmable one by adding two 
registers specifying the values of the constants a and b. 

Let us now address the issue of how to find the param- 
eters a and b entering in © mapping line L and L'. For 
this, we consider the lines defined as 

L : n\q + n 2 p — mod 2N 
L' : q +p' = n 3 mod 2N. 

For simplicity, we consider the case where at least one 
of the parameters rii is an odd number (the other case 




FIG. 2: Programmable gate array to evaluate the average 
Wigner function over a tilted line in phase space. The cat map 
is parameterized by two integers, which can be programmed 
in auxiliary registers (so the array can be made fully pro- 
grammable) . 

can be treated similarly). It is worth mentioning that 
the program state for L' can be efficiently prepared. The 
mapping between the two lines is accomplished by using 
a cat map as in JSJ) with the parameters given by 

n 2 a = 1 — m mod 2N 

b = l + n 2 mod 2N. (9) 

This method can be generalized to evaluate the aver- 
age value of the Wigner function over tilted rectangular 
regions. For this purpose, we can construct the program 
state for a simple rectangular region (defined by the con- 
ditions: qi < q < (72, Pi < P < P2)- Then, we can 
map this region into a tilted region by using the strategy 
described above. 



III. PROGRAMMABLE TOMOGRAPHY OF 
THE KIRKWOOD DISTRIBUTION 

The Kirkwood function is a phase space distribution 
whose use is probably less common. It was first pro- 
posed by Kirkwood 14], and used in quantum statistics. 
It displays different phase space features of a quantum 
state in phase space and has the advantage of being di- 
rectly linked to the matrix elements of the density matrix 
in a mixed representation. On the other hand the Kirk- 
wood function is a complex number even for hermitian 
operators. Having defined the unitary basis of operators 
{$1^5 ?>P = 0, .., N — 1} in a Hilbert space of dim A, 
the discrete Kirkwood function of a density operator p is 
defined as: 

where \q) and \p) are position and momentum eigen- 
states, respectively. The circuit that implements the 
measurement of the Kirkwood distribution is also based 
on the scattering circuit. It can be seen in Figure |3 
where the swap gate acts as 5[|^)i|$) 2 ] = |$)i| v I / )2- It 
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FIG. 3: Programmable gate array to measure the Kirkwood 
function. The first register is the ancillary (probe) qubit, the 
second and the third registers are the program states. The 
FT gate is the quantum Fourier transform, and the swap gate 
acts as S[|*)i|$) 2 ] = [$)i|*)a. 



is simple to show that if the program state (second and 
third registers) are computational states specifying the 
position and momentum coordinates, this circuit allows 
to measure the Kirkwood distribution at any point of the 
phase space. By measuring the expectation values of (a z ) 
and (a y ) for the ancillary qubit we obtain the Kirkwood 



distribution as 



Tr \p)(p\p\q)(q\ = K(q,p) 



(11) 



IV. PROGRAMMABLE TOMOGRAPHY OF 
THE HUSIMI DISTRIBUTION 

The Husimi function is a well known alternative dis- 
tribution in phase space, which is based on the use of 
minimum uncertainty wave packets \a). The Husimi dis- 
tribution is the expectation value of the density matrix in 
the coherent state \a). This is a positive quantity that, in 
the continuous case, graphically displays the phase space 
contents of the state in a region of area h. In the dis- 
crete case the coherent states can also be defined 0] 
(we provide below an efficient scheme for their prepara- 
tion). These wave packets define an over-complete basis 
{|a),ev = (q,p),q,p — 0,..N— I}. In terms of this basis 
the Husimi distribution is defined as 



H(a) = ^Tr 



\a)(a\p 



(12) 



The programmable circuit that implements the mea- 
surement of the Husimi distribution is also based on the 
scattering circuit (in particular, it is a straightforward 
application of the ideas proposed in |l7j and used in the 
experiment presented in [l8lp. In Figure 0] we can see a 
representation of the algorithm, which uses an ancillary 
(probe) qubit, and a program register prepared in the 
state | a) (a coherent state centered at the point where 
we want to evaluate the Husimi distribution). It is easy 
to show that the circuit is such that 



Tr \a)(a\p = NH(a). 



(13) 



The algorithm would only be useful if an efficient method 
can be devised to prepare the state \a). We devote the 
next section to a complete description of this subroutine. 
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FIG. 4: Programmable circuit to measure the Husimi distri- 
bution function. The first register is the ancillary (probe) 
qubit, the second register is the program state (a coherent 
sate centered at the point of interest in phase space) and the 
third one is the system of interest. 



V. EFFICIENT ALGORITHM FOR THE 
GENERATION OF COHERENT STATES 

A. Discrete coherent states 

In the continuous case, coherent states can be defined 
as phase space translations acting on the ground state 
of the harmonic oscillator Hamiltonian. In the discrete 
case with periodic boundary conditions, the harmonic 
oscillator can be replaced by the Harper Hamiltonian: 

u + W v + v^ 



|$o> =£o|*o), 



(14) 
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which ensures the proper periodicity conditions. Coher- 
ent states can now be defined |16| as discrete translations 
on |$ ) 



(15) 



where T(a) = U q V p e 2 ^' pq are the phase space translation 
operators. An alternative definition, yielding a continu- 
ous distribution with analytic properties |l9| is given 
by: 



e -*TMW-*-»0|n). (16) 

These states are almost indistinguishable from (|15f) as N 
grows, and both are periodic wave packets occupying a 
minimum uncertainty area 1/N in phase space. Their de- 
tailed structure, showing the extremely small differences 
at small values is shown in Figure [5] 

The classical analogue of the Harper Hamiltonian i|14|) 
is H = ^(sin 2 7rQ-|-sin 2 ttP). This gives rise to the follow- 
ing classical map equations (for a small time step 7/2-71") 

Q' = Q-7sin(27rP) modi (17) 
P' = P + 7sin(27rQ') modi. (18) 
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FIG. 5: Population (in the coordinate basis) of both the con- 
tinuous 1161 and the discrete 1141 coherent states, centered 
at the origin of the phase space. The difference between both 
definitions tends to zero in the large N limit (here, N = 64). 



FIG. 6: Population (in the position basis) of the eigenstates 
of the Harper Hamiltonian and the Harper kicked map, lA{"i), 
as a function of 7. We can appreciate that for the smallest 
values of 7 these states are almost identical (N = 64). 



For infinitesimal 7 we obtain Hamilton equations, and the 
conservation of energy leads to integrable behavior. Our 
strategy will be to quantize the map equation for small 7 
as a way to obtain a unitary operator with an eigenstate 
very close to |$o). The map belongs to the well known 
family of the kicked maps and the unitary operator cor- 
responding to its quantization is obtained as a product 
of two operators representing a potential and a kinetic 
kick. These two operators are respectively diagonal in 
position and momentum basis. They can be efficiently 
implemented by means of a quantum network consisting 
of two controlled phases interposed by the Fourier trans- 
form: 



W( 7 ) = M[FT]^K[FT], 



(19) 



where FT is the iV-dimensional Fourier transform and 
the operators M y K represent the potential and kinetic 
kicks respectively: 



M\q) = e -^ Ncos ^\q) 



K\p) 



\P)- 



(20) 
(21) 



The crucial feature of this unitary operator is that as 
7 <C 1 its eigenstates become those of Harper Hamilto- 
nian. Hence, U (7) has a coherent state as one of its eigen- 
states for small values of 7. In Figure we show how the 
eigenstate of the map lfr3|l converges to the ground state 
|$o) as 7 -> 0. 

From the above discussion is clear that what we need is 
a method to efficiently prepare an eigenstate of the uni- 
tary operator (|19fl . For this, we will use the well known 
phase estimation algorithm [2£j, l2ll 12^ to filter an ini- 
tial state which is approximately localized near the phase 
space origin. 



B. Algorithm for phase estimation 

As the phase estimation algorithm [23, 12H 12^ , is an 
essential part of our construction we briefly review its 
operation. The circuit is reproduced in Figure Its 
operation in the phase estimation mode requires that an 
eigenstate of W( 7 ) be supplied to the lower register. Then 
a measurement performed in the upper registers yields 
a rational approximation to the eigenphase, which im- 
proves as the size of the upper register increases. Here 
we are more interested in the use of the circuit as a fil- 
ter in which case a state approximating an eigenstate is 
fed to the lower register. If this state is expanded as 
1^) = J2 a c a\ u o') the measurement in the upper register 
yields a distribution of phases with probabilities propor- 
tional to |c Q | 2 . Furthermore, if the number of qubits in 
the first register is such that the value of the phase can 
be exactly determined, then the final state of the second 
register is the eigenstate corresponding to that phase. 
Thus, it is clear that if the initial state \ty) is near \&o), 
the application of this circuit would provide us with the 
desired state with high probability (see below). 
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FIG. 7: The circuit for the phase estimation algorithm. When 
the state of the lower register is an eigenstate of U, the mea- 
surement of the upper register reveals an t-bit approximation 
to the corresponding eigenphase. 

However, the typical situation is when the number of 
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qubits in the first register only allows us to make an ap- 
proximation to the real value of the phase. Thus, if the 
system has n qubits and the number of qubits in the first 
register is 

t = n + log 2 [2 + 1] , (22) 

the phase estimation algorithm gives an approximation 
to the phase corresponding to the eigenstate \u), with 
probability bounded by 

Pu = \c u \ 2 (l-e). (23) 

This clearly differs from the ideal case, since now we do 
not obtain the exact value of the phase. Therefore the 
above probability corresponds to the measurement of an 
integer k such that 4 is the best i-bit estimate to ip (with 
an error bounded by |2*~™ — 1|). After the measurement, 
the state of the second register will be a linear combi- 
nation of all the eigenstates of U that is close to the 
corresponding eigenstate (see below). 
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FIG. 8: The Wigner function of a square state. The state 
is located in position and partially localized in momentum 
(N=64). Horizontal (vertical) axis corresponds to position 
(momentum) basis. The color convention is such that positive 
(negative) values of the Wigner function correspond to black 
(white) regions. 



C. Algorithm for the generation of coherent states 

As mentioned above, the algorithm for the prepara- 
tion of coherent states consists of the application of 
the phase estimation algorithm to filter an initial state 
(which should be itself a well localized state near the ori- 
gin of phase space). We will now describe in detail the 
two necessary ingredients for the efficient implementation 
of the algorithm: i) the preparation of the initial state 
and, ii) the efficient implementation of the ctvl-W (7) 
gates, required for the phase estimation to be applica- 
ble. 



1. Initial state preparation 

This is indeed the simplest part of the algorithm. In 
fact, an easily preparable candidate for the initial state 
is what we could denote as a "square state", defined as 
an equally weighted superposition of the first v^/V states 
of the computational basis: 

w—l 

i*°> = ^Eif -9} modiV ' (24) 

q=0 

where N is the Hilbert space dimension and w — [\ZiVj. 
This state is strictly localized in position in a region of 
width w around the origin but, because of diffraction, is 
only partially localized in momentum. This can be seen 
by analyzing its Wigner representation shown in Figure|Hl 
The important features of this state are that it has strong 
overlap with a coherent state localized at the origin of 
phase space (which is our target state). In the limit of 
large N, the overlap tends to a value of 0.94. Also, it 
is simple to show that the square state can be efficiently 



prepared: Starting from |0)® n we simply need to apply 
Hadamard gates to the n/2 least significant qubits we 
obtain a state centered at the phase space point (w/2, 0). 
Centering the state at the origin requires a shift, which 
can be implemented efficiently. 

2. Efficient implementation of the ctri-W (7) gates 

The phase estimation algorithm requires that the pow- 
ers of the operator U (7) be implemented efficiently. This 
is in general not the case. However, we can get around 
this problem by using the fact that the Harper map has 
a well defined semi-classical limit. This allows us to per- 
form its iteration by using a reliable semi-classical ap- 
proximation. Thus, for 7 <C 1 (and for N large) the 
powers needed can be approximated as follows: 

« 2 *( 7 )~W(2*7). (25) 

This approximation also relies on the assumption that 
the initial state is localized in a region of the phase space 
where the map is regular (this is satisfied by the square 
state). If this approximation is valid, then the powers of 
the unitary operator can be implemented efficiently us- 
ing the same quantum networks required to implement 
the operator itself (to implement a power of U we simply 
use a different parameter 7). This approach has a clear 
limitation: Each power increases the value of 7 and for 
large values of 7 the map ceases to be integrable. In such 
case, its spectral properties become very different from 
those of the Harper Hamiltonian. Our goal then is to 
propagate the map for a time long enough to resolve the 
ground state from its neighbors without violating the ap- 
proximation l|25|) . These two issues (evolving accurately 
and resolving the spectrum) should be studied jointly. 
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FIG. 9: The fidelity (overlap between the state evolved with 
the exact and approximated unitary operators) as a function 
of the number of qubits in the first register. 2 ts defines the 
maximum time for which 12-") I is valid. 



FIG. 10: Dependence of the phase difference Aip with 7. A 
linear dependence is found for values of the Hilbert space 
dimension N such that 7V7 < 0.6. 



But to make our presentation clear we can first analyze 
them separately. 

To examine the accuracy of the approximate evolu- 
tion we can study the fidelity, defined as the absolute 
value of the overlap between the square state propa- 
gated with the exact and approximate evolution (i.e. 
F = ^\(iJj eX act\?papprox)\) ■ In Figure El we plot this as 
a function of the number of qubits in the first register. 
We find that the fidelity remains close to unity up to a 
sharply defined time T s (T s = 2 ts ), after which it drops 
abruptly. This time, T s , defines the allowed number of it- 
erations of the unitary compatible with a given accuracy. 
Therefore, it sets a bound to the maximum number of 
qubits we can include in the first register of the phase es- 
timation algorithm for the semi-classical approximation 
to remain valid. 

To analyze the resolution required to resolve the spec- 
trum of the unitary operator we should consider the mini- 
mum difference between neighboring eigenphases of U(^) 
(denoted as Atp) . This quantity determines the minimum 
number of qubits in the first register, needed to resolve 
the ground state. Thus, for this purpose, we would need 
a number of qubits n which should be at least equal to 

„ = log a (J-). (26) 
Aip 

An important question is how does this number scales 
with the dimension of the Hilbert space of the system. 
This can be determined by analyzing the dependence of 
Aip with 7 and N. This is done in Figure ITUI where we 
show that the phase difference has a linear dependence 
with 7, at least for values of N and 7 such that A^7 < 0.6. 
Therefore, the number of qubits required to resolve the 
spectrum will scale logarithmically with the dimension- 
ality of the system's Hilbert space. 

As mentioned above, the behavior of the fidelity as a 



function of the parameter 7 should be analyzed jointly 
with the minimum number of qubits needed in the first 
register, to achieve the required spectral resolution ac- 
cording to eq. (|26p. In Figure ITT1 we display the curves 
coming from each of these requirements. To achieve high 
enough fidelity the value of 7 and of the number of qubits 
in the first register must be below the lower curve. On 
the one hand to achieve the required spectral resolution 
one needs the value of 7 and the number of qubits to 
lie above the second curve. This seem to be a problem, 
but can be easily solved taking into account the follow- 
ing observations. Thus, we notice that: i) the two lines 
are approximately parallel (for all values of N, as long as 
Nj < 0.6) and ii) the two parallel lines are simply shifted 
away from each other by about three qubits. Hence, the 
semi-classical approximation can be used up to a power 
of W(7) given by M("/) 2 ™ with n — log 2 (^). To iterate 
the map further, as required to achieve enough spectral 
resolution we can apply the remaining powers of U (7) 
as products of its precedents (which were implemented 
using the semi-classical approximation). As the distance 
between curves is fixed, we would only need to use this 
trick a fixed (A-independent) number of times. Doing 
this, it is possible to achieve the required accuracy and 
spectral resolution simultaneously. 

It is worth mentioning that the precision needed to re- 
solve the ground state is not the only condition that im- 
poses a lower bound on the number of qubits of the first 
register. Thus, in principle if we want to get the desired 
coherent state with a reasonable probability we need the 
number of qubits to obey the relation fixed by equation 
(|23|l . Suppose that we impose a value of e = 1/4, which 
corresponds to a probability for preparing the right co- 
herent state of about po = 0.70 (this is computed taking 
into account that the initial square state is such that 
I Co 1 2 =0.94). For this value of e, equation (|23fl implies 
that the first register should have at least two more qubits 
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In Figure ITU we show how the approximation can be im- 
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FIG. 11: Dependence of the parameter 7 as a function of 
the number of qubits in the first register. The dashed line 
represents the condition of high fidelity in the evolution (the 
allowed parameter region lies below this curve). The solid 
line corresponds to the condition of good spectral resolution 
(the allowed parameter region lies above that curve) . The two 
lines are shifted from each other by about three qubits. 



than the system's register. This is a lower bound for the 
dimension of the first register. 



3. Final remarks on the preparation of coherent states 




FIG. 12: Wigner function for four of the quantum states of the 
system's register that are generated by the phase estimation 
algorithm when the detected value of the first register falls 
in the peak of the probability distribution. Horizontal (verti- 
cal) axis corresponds to position (momentum) basis. Labeling 
these states horn left to right and from top to bottom as a), 
b), c) and d), we can see in Figure [TBI their representation in 
the computational basis. 



In summary, the algorithm to prepare a coherent 
state consists of the following steps: i) preparation of 
a "square" state, ii) selection of the parameter 7 for the 
evolution operator and the corresponding determination 
of the number of qubits to be used in the first register 
of the phase estimation algorithm, iii) run the phase es- 
timation algorithm efficiently implementing the powers 
of the operator U{p/) in an approximate way, iv) from 
the peaked distribution of results for the phase, we dis- 
cover the one associated with the coherent state, when 
this phase is detected the desired coherent state has been 
prepared in the system's register, v) after obtaining a co- 
herent state centered at the origin, one can translate it to 
any point of the phase space using phase space displace- 
ment operator, which can be efficiently implemented as 
in the circuit of Figure ^ 

To see the algorithm in action we performed a few nu- 
merical simulations. In Figures |12I we show the Wigner 
function of four of the quantum states of the system's 
register that fall under the peak of the probability distri- 
bution for the first register (we used N — 64 and showed 
in Figure 1131 the probability distribution for the same 
states in the computational basis). It is clear that any of 
such states is a good approximation to a coherent state 
(the fact that this stage of the algorithm works as a fil- 
ter can be appreciated by comparing the initial square 
state shown in Figure [S] and those shown in Figure IT5f > . 



proved even further by various means. In fact, we could 
include several stages of filtering each one of which would 
considerably improve the quality of the final state. 




FIG. 13: Probability distribution in the computational basis 
for four of the quantum states of the system's register that 
are prepared when the detected value of the first register falls 
in the peak of the distribution. The Wigner function of such 
states is shown in Figure [T21 
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FIG. 14: Probability distribution in the computational basis 
for the states produced after different number of iterations of 
the filtering algorithm. Each iteration improves the quality 
of the state, which becomes closer to a true coherent state. 



of the circuits is based on the fact that operations such 
as phase space translations, reflections and the Fourier 
transform are efficiently implementable. For the case 
of Wigner and Kirkwood distributions, the efficiency is 
solely based on this fact (however, it is worth pointing 
out that, contrary to what happens with the Wigner and 
Husimi distributions, the evaluation of a typical value of 
the Kirkwood distribution of a pure state would require 
exponential precision due to the factor of TV absent in 1|11|) 
as compared with Qj). The evaluation of the Husimi 
distribution requires the use of a subroutine preparing 
coherent states. We presented a method achieving this 
goal, which consists of a variation of the phase estimation 
algorithm with an appropriately chosen initial state. In 
this case the efficiency requires not only a good guess for 
the initial state (which is indeed easily done) but also the 
possibility of efficiently implementing powers of the uni- 
tary map whose eigenstate is close to a coherent state. 
In our case, this can be done by using a semi-classical 
approximation for this operator. 
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